#ifndef I_LIB_IP_I_INTEGRAL_H
#define I_LIB_IP_I_INTEGRAL_H

#include "../core/i_blas.h"

namespace idl
{
    inline void i_integral_image_32f_inplace(float *const *image, int width, int height)
    {
        int x, y;
        float acc;
        //first pixel
        image[0][0] = image[0][0];
        //first row:
        for (x = 1; x < width; ++x)
        {
            image[0][x] = (image[0][x] + image[0][x - 1]);
        }
        //handle remaining rows:
        for (y = 1; y < height; ++y)
        {
            acc = image[y][0];
            image[y][0] = image[y - 1][0] + acc;
            for (x = 1; x < width; ++x)
            {
                acc += image[y][x];
                image[y][x] = image[y - 1][x] + acc;
            }
        }
    }

    inline void i_integral_image_32f(const float *const *image, 
        float **integral, 
        int width, 
        int height)
    {
        int x, y;
        float acc;
        //first pixel
        integral[0][0] = image[0][0];
        //first row:
        for (x = 1; x < width; ++x)
        {
            integral[0][x] = (image[0][x] + integral[0][x - 1]);
        }
        //handle remaining rows:
        for (y = 1; y < height; ++y)
        {
            acc = image[y][0];
            integral[y][0] = integral[y - 1][0] + acc;
            for (x = 1; x < width; ++x)
            {
                acc += image[y][x];
                integral[y][x] = integral[y - 1][x] + acc;
            }
        }
    }
  
    inline void i_integral_image_32i(const int *const *image, int **integral, int width, int height)
    {
        int x, y;
        //first row:
        int acc = 0, *ptri = NULL;
        const int *cptrs = NULL;
        for (x = 0; x < width; ++x)
        {
            acc += image[0][x];
            integral[0][x] = acc;
        }
        //remaining rows:
        for (y = 1; y < height; ++y)
        {
            acc = 0;
            cptrs = image[y];
            ptri = integral[y - 1];
            for (x = 0; x < width; ++x)
            {
                acc += cptrs[x];
                integral[y][x] = acc + ptri[x];
            }
        }
    }

    template <typename S, typename T>
    inline void i_integral_image(const S *const *image, T **integral, int width, int height)
    {
        int x, y;
        //first row:
        T acc = (T)0.0, *ptri = NULL;
        const S *cptrs = NULL;
        for (x = 0; x < width; ++x)
        {
            acc += (T)image[0][x];
            integral[0][x] = acc;
        }
        //remaining rows:
        for (y = 1; y < height; ++y)
        {
            acc = (T)0.0;
            cptrs = image[y];
            ptri  = integral[y-1];
            for (x = 0; x < width; ++x)
            {
                acc += (T)cptrs[x];
                integral[y][x] = acc + ptri[x];
            }
        }
    }

    /*Computes the sum of pixels within the rectangle specified by the upper-left corner (x, y)*/
    template <typename T>
    inline T i_integral_rect(const T *const *integral, 
        int x, int y, 
        int rect_width, int rect_height, 
        int integral_width, int integral_height)
    {
        int x_1 = i_min(i_max(0, x - 1), integral_width-1);
        int y_1 = i_min(i_max(0, y - 1), integral_height-1);
        int x_2 = i_max(1, i_min(x + rect_width, integral_width)) - 1;
        int y_2 = i_max(1, i_min(y + rect_height, integral_height)) - 1;
        T a = integral[y_2][x_2];
        T b = integral[y_2][x_1];
        T c = integral[y_1][x_2];
        T d = integral[y_1][x_1];
        return i_max(a - b - c + d, (T)0.0);
    }

    /*Computes the sum of pixels within the rectangle specified by the upper-left corner (x, y), the actual area (# of pixels) 
      of the valid rectangle is returned. With boundary conditions, the actual number of pixels may not equal to rect_width x rect_height*/
    template <typename T>
    inline T i_integral_rect(const T *const *integral, 
        int x, int y, 
        int rect_width, int rect_height, 
        int integral_width, int integral_height, 
        int& valid_rect_area)
    {
        int x_1 = i_min(i_max(0, x - 1), integral_width - 1);
        int y_1 = i_min(i_max(0, y - 1), integral_height - 1);
        int x_2 = i_max(1, i_min(x + rect_width, integral_width)) - 1;
        int y_2 = i_max(1, i_min(y + rect_height, integral_height)) - 1;

        /*the actual area (# of pixels) of the valid rectangle region*/
        valid_rect_area = i_max((y_2 - y_1) * (x_2 - x_1), 1); 
        
        T a = integral[y_2][x_2];
        T b = integral[y_2][x_1];
        T c = integral[y_1][x_2];
        T d = integral[y_1][x_1];
        return i_max(a - b - c + d, (T)0.0);
    }

    /*Computes the sum of pixels within the rectangle specified by the upper-left corner (x, y), this routine assumes the input x, y, rect_width, rect_height
      have been validated therefore skips the step of checking boundary conditions.*/
    template <typename T>
    inline T i_integral_rect_safe(const T *const *integral, 
        int x, int y, 
        int rect_width, int rect_height, 
        int integral_width, int integral_height)
    {
        int x_1 = x - 1;  
        int y_1 = y - 1;  
        int x_2 = x + rect_width - 1; 
        int y_2 = y + rect_height - 1;
        T a = integral[y_2][x_2];
        T b = integral[y_2][x_1];
        T c = integral[y_1][x_2];
        T d = integral[y_1][x_1];
        return i_max(a - b - c + d, (T)0);
    }

} /* namespace idl */

#endif